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I.  Introduction 

■w 

We  present  a  numerical  study  of  the  blow-up  of  u,  =  +  u''.  I  his  is  one  of 

a  large  class  of  nonlinear  evolution  equations  with  scale-invariant  structure  and 
blowing-up  solutions.  Other  examples  include  reaction  diffusion  equations  such 
as  M,  -  Au  =  u''  or  m,  -  Au  =  e“.  which  arise  in  models  of  combustion 
(e.g.  [21],  [22]).  and  the  nonlinear  Schrodinger  equation  iu,  -  Am  =  |m|''’  'm. 
which  ari.ses  in  plasma  physics  and  nonlinear  optics  (e.g.  |24|.  |25|).  f  he  blowing- 
up  solutions  of  such  equations  have  in  common  the  properties  that  (i)  the 
singularities  are  isolated,  and  (ii)  the  singularities  have  a  characteristic  structure, 
which  may  or  may  not  be  directly  linked  to  the  scaling  properties  of  the  equation. 
For  .solutions  that  do  develop  singularities  it  is  often  of  interest  to  study  the  local 
features  of  the  blow-up.  Such  a  study  may  be  useful  not  only  for  direct 
comparison  to  the  phenomenon  being  modelled,  but  also  for  extending  the 
solution  beyond  the  singular  time,  or  for  matching  it  to  the  solution  of  a  different 
equation  which  applies  near  the  singularity. 

Several  authors  have  attempted  to  calculate  the  local  character  of  the  blow-up 
of  M,  =  M,,  +  u’’  numerically  (e.g.  [6].  [13],  [20]).  This  is  a  sensitive  problem, 
since  most  methods  for  .solving  evolution  equations  lose  accuracv  as  the  solution 
becomes  large.  Two  novel  approaches  have  recently  been  introduced,  in  some¬ 
what  different  contexts;  Chorin  used  an  algorithm  based  on  rescaling  and  mesh 
refinement  to  study  the  three-dimensional  Euler  and  Navier-Stokes  equations  in 
[7j;  and  a  method  based  on  continuous-in-time  rescaling  has  been  applied  by 
LeMesurier.  McLaughlin,  Papanicolaou.  P.-L.  Sulem.  and  C'.  Suleni  to  study  the 
nonlinear  Schrodinger  equation  in  [24],  125]. 

Our  approach  differs  from  those  just  cited  in  the  following  was.  L’pon 
re.scaling  to  resolve  the  appearing  singularity.  C'horin's  method  concentrates  on 
an  increasingly  small  physical  domain,  enforcing  periodic  boundary  <'onditions 
closer  and  closer  (in  unsealed  distance)  to  the  singularity.  The  method  of 
LeMesurier  et  al.  uses  a  fixed  mesh  in  physical  space  which  is  spread  apart  by 
rescaling,  so  that  accuracy  is  inevitably  lost  far  from  the  singularity.  In  contrast, 
using  our  me.sh  refinement  strategy  we  are  able  to  compute  accurately  over  the 
entire  physical  interval  even  as  the  solution  grows  in  magnitude  from  ()( I )  to 
0(10*').  The  main  idea  is  this:  we  step  the  solution  forward  until  its  maximum 
value  reaches  a  preset  threshold.  Where  the  resulting  function  is  large  the 
solution  is  rescaled  to  make  it  small  again.  Since  scaling  stretches  the  spatial 
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variable,  extra  grid  points  are  added  to  maintain  accuracy.  The  rescaled  solution 
is  then  stepped  forward  until  its  maximum  value  reaches  the  threshold  value,  at 
which  juncture  a  further  rescaling  takes  place,  etc.  In  effect,  our  procedure  solves 
the  equation  with  a  varying  mesh  width  and  time  step  that  are  linked,  at  each 
point  in  space-time,  to  the  magnitude  of  the  solution.  Though  very  similar  to 
Chorin’s  method  in  spirit,  ours  has  the  advantage  that  the  boundary  conditions 
for  the  rescaled  problems  are  handled  in  a  manner  that  is  consistent  with  the 
underlying  evolution  equation.  Our  approach  to  mesh  refinement  and  multiple 
grids  is  similar  to  one  which  has  been  used  for  solving  first-order  hyperbolic 
systems  in  one  or  more  space  dimensions;  see  e.g.  [4].  In  fact,  we  implemented 
the  algorithm  by  modifying  a  code  originally  developed  to  solve  hyperbolic 
systems  in  one  space  dimension. 

Though  the  method  is  obviously  more  general,  we  have  applied  it  only  to  the 
semilinear  heat  equation 

(1.1)  u,  =  M,,  -  u’’.  P  >  1. 

on  the  interval  -1  <  .v  <  1.  with  a  Dirichlet  boundary  condition  u{-\.i)  = 
u(\.  t)  =  0.  Attention  is  further  restricted  to  initial  data  <)(.\ )  such  that 

(1.2)  <f>  >  0.  (p(x)  =  <{)(-.v).  <  0  for  .V  ^  0. 

for  which  the  solution  of  (1.1)  is  positive,  symmetric,  and  radially  decreasing.  A 
lot  is  known  about  how  solutions  of  this  equation  blow  up;  .see  e.g.  (1],  [2],  [5]. 
[10]-(19],  [23],  [26].  In  particular,  one  knows  that 

(1.3)  lim(7'- ‘’M(|/r^./)  =  (/7  -  1)  '  ''■  " 
r  f  r 


uniformly  for  |4|  <  C,  where  T  is  the  blow-up  time;  see  [14],  [16].  This  gives  the 
behavior  in  any  space-time  parabola  |.y|’  <  C{T  -  t)  based  at  the  blow-up  point. 
It  is  natural  to  ask  what  happens  beyond  these  parabolas.  For  example,  what  is 
the  asymptotic  shape  of  the  curve  where  {T  -  "w  is  constant? 

In  [13],  [14],  Galaktionov  and  Po.sashkov  u.se  a  formal  argument  adapted  from 
[20]  to  derive  the  ansatz 


(1.4) 


I 

u(x.t)  ~{T- l)’’  ' 


(p-l)  + 


ip  -  1)' 
4/7 


(T 


z)|iog(r  -  /)| 


This  is  consistent  with  (1.3).  since  the  .second  term  in  the  bracket  tends  to  zero  as 
t  T  with  I  =  x/  -  t  fixed.  It  suggests  that  the  curves  (T  -  "»  =  y 

are  a.symptotically  of  the  form  =  o(y )(7  -  z)|log(  T  -  /  )|.  (An  analogous 
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conjecture  concerning  the  blow-up  of  m,  -  is  presented  in  [9].)  Some 

preliminary  analytical  results  tending  to  confirm  (1.4)  are  given  in  [10]  and  [14], 
but  they  fall  far  short  of  a  full  proof.  Our  numerical  calculations  give  sufficient 
detail  of  the  behavior  near  blow-up  to  allow  us  to  test  (1.4).  The  calculated 
solution  is  in  fact  in  excellent  agreement  with  the  formula,  leaving  little  doubt  in 
our  minds  about  the  validity  of  this  conjecture. 


2.  The  Algorithm 

The  basis  of  our  algorithm  is  the  following  scale  invariance  of  the  equation 

(1.1):  if  u(x.  I)  solves  it.  then  so  does 


(2.1) 


Uyiy.  t)  =  y’'''''  yV) 


for  any  y  >  0-  By  choosing  y  to  be  small  when  u  is  large,  one  can  keep  the 
rescaled  solution  m,  bounded.  Thus  it  is  easier  to  solve  for  than  it  is  to 
compute  u  directly.  Note  however  that  both  space  and  time  are  stretched  by  the 
scaling;  if  u  is  defined  for  -1  <  x  <  1  and  0  <  /  <  7.  then  the  domain  of  is 
-y  ‘  <>'<y“'.0<T<Y‘"7.  This  is  the  price  one  pays  for  the  advantages  of 
rescaling.  Computationally,  if  u  is  defined  on  a  grid  of  mesh  width  A.v.  then  (2.1 ) 
defines  only  on  a  grid  of  mesh  width  y "  ’  A.v.  A  loss  of  accuracy  is  avoided  by 
introducing  additional  points  to  the  grid  on  which  is  defined.  Our  algorithm 
maintains  both  the  original  w  and  the  rescaled  u^.  each  defined  on  a  separate 
grid,  and  steps  each  forward  in  a  time-accurate  way.  Since  the  scaling  (2.1) 
stretches  time  as  well  as  space,  most  of  the  computational  effort  goes  into  the 
advancing  of  Uy.  Actually,  our  algorithm  has  an  iterative  structure,  so  that  at  the 
/c-th  iteration  we  are  maintaining  not  just  one  rescaled  solution  but  k  of  them, 
corresponding  to  y  =  X,  A*,  where  X  is  a  fixed  .scaling  parameter.  To 

avoid  unnecessary  computation  we  rescale  only  where  u  is  large,  in  such  a  way 
that  the  finest  rescaled  solution  (m.^  with  y  =  X*  in  the  A  iteration)  stays  bounded 
away  from  0  as  well  as  oc. 

To  determine  the  algorithm  one  must  fix  three  parameters: 

X  =  scale  factor. 


(2.2)  M  =  maximum  height  before  re.scaling. 

a  =  parameter  controlling  width  of  the  interval  to  be  rescaled 
(which  is  where  aM  £  u  ^  M). 


They  should  be  chosen  so  that 

(2.3)  X  ‘  >  1  is  a  small  integer,  and  0  <  a  <  1 . 


Typical  choices  are  X  = 


M  =  2/2 .  0=1.  The  algorithm  computes  succes- 
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sively  a  sequence  of  functions  where 

=  ^-th  rescaled  solution, 

(2.4)  }\  =  A-th  rescaled  spatial  variable. 

=  clock  measuring  (rescaled)  time  for 

The  time  (according  to  the  clock  t^)  at  which  is  re.scaled  to  yield  ,  will  be 
denoted  by  t^*.  and  the  interval  which  is  re.scaled  will  be  (  \\  .  \\  ).  All  these 
quantities  will  be  defined  more  precisely  below.  The  initial  index  A  =  0  corre¬ 
sponds  to  the  “real”  solution  u  as  a  function  of  “real"  space  .v  and  "real"  time  i: 

(2.5)  «(,=  «.  ,b)  =  V.  T„  =  t. 

The  initial  phase  of  the  algorithm  .simply  integrates  the  equation  (1.1)  until 
the  maximum  amplitude  reaches  M.  (We  assume  that  the  initial  data  satisfy 
<t>  ^  M,  and  that  the  corresponding  solution  u  does  indeed  blow  up.  This  is  true, 
for  example,  if  <))  =  c(l  +  cos(w.v))  with  c  sufficiently  large,  and  A/  >  2(.)  We 
use  the  forward  Euler  finite  difVerence  scheme 


(2.6) 


»(  v,.  i„ .  ,)  =  a(  v,.  !„) 
At 


— —  jt/(.v,  |.  /„)-2tt(.v,.  /„)  -  »(  A,. ,.  .'„)]  +  A/  •  a '(a,.  r„  I 
Avr* 


which  is  first-order  accurate  in  time  and  second-order  in  space.  Typical  choices 
are  Ax  =  .005  and  At  =  i(A.v)’.  In  this  initial  phase  of  the  algorithm  the 
Dirichlet  condition  «  =  0  is  used  at  the  endpoints  .\  =  ±1.  The  solution  is 
integrated  until  the  first  time  step  when  ||m( •.  t„)|l.^  >  Af.  Then  it  is  linearly 
interpolated  in  time,  using  two  time  levels,  to  obtain  a  time  t,,*  with  !„  -  At  s 
''i*  ^  T,  such  that 

(2.7)  ||«(-.r„*)t  =  .W. 

The  interval  to  be  rescaled  at  the  next  stage.  (  f,, .  i,,' ).  is  essentially  the  .set  w  here 
«(*.  V)  ^  aM.  It  is  convenient  however  to  let  i„'  he  grid  points,  so  they  are 
defined  by 


(2.8) 


m(  h,  -  A.v.  T„*  )  <  ttA/  ^  m(  r„  ,  T,,*  ). 
“(.'»  •  V  )  ^  «.A/  >  w(  i',;  +  Aa.  t,*  ). 


Of  course,  perfect  arithmetic  would  yield  v„  =  -  We  do  not  enforce  tins 
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Figure  1.  The  graph  of  w„  jusi  pnor  to  rescaling. 


condition,  so  that  symmetry  may  be  used  as  an  indicator  of  the  accuracy  of  the 
computation.  Figure  1  depicts  the  graph  of  u  at  the  end  of  this  initial  phase,  just 
prior  to  rescaling. 

The  first  “rescaled  solution”  u,  is  related  to  u  by 

(2.9)  M,(  i-,.  T, )  =  A’'*'’  "m(A  V|.  T|*  +  A’t,  ). 

We  want  u  to  be  evaluated  on  the  right  of  (2.9)  only  where  u  s  aM:  therefore  \ y 
is  restricted  to  the  interval 

(2.10)  A  ^  ‘i„  g  !  [  ^  A  '  V|,'  . 

The  maximum  value  of  u^  at  its  initial  time  t,  =  0  (corresponding  to  i  =  t,,*  )  is 
reduced  from  M  =  ||m(',  to 

I1mi(-.0)||^  =  A-"'^  "M  <  M: 

this  is  the  purpose  of  rescaling.  Due  to  the  scale  invariance  property  (2.1).  m, 
solves  the  same  equation  as  u.  with  respect  to  its  (rescaled)  “space"  and  “time" 
variables  y,  =  a/A,  t,  =  (/  —  t,*)/A‘: 

JL  11  -  r 


(2.11) 
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This  is  what  lies  at  the  foundation  of  our  algorithm;  the  difference  scheme  (2.6) 
originally  introduced  for  stepping  u  =  forward  in  time  can  also  be  used  to 
solve  for  m,  (and,  eventually,  for  each  successive  ). 

The  computation  of  u^  requires  initial  data  wd'-O)  boundary  data 
u^(\  ‘v'o^ ,  Tj).  The  former  are  obtained  by  rescaling  u: 

(2.12a)  =  A'"*'’  V). 

The  latter  are  obtained  using  the  coarse  mesh  appro.ximation  method  (see  [8]): 
the  boundary  condition  for  the  refined  problem  is  determined  by  the  solution 
previously  computed  in  the  same  region  using  the  coarser  mesh.  Specifically,  the 
right  side  of  equation 

(2.12b)  u,(X-V(r-  VAt,)  =  A-/'"  “«(>•»■  ,  T,*  +  A^^At,) 

is  obtained  by  applying  the  forward  Euler  difference  scheme  to  u  with  a  time  step 
A-  At,  (which  is  smaller  than  the  full  time  step  At^  =  A/  used  to  advance  u  on  the 
coarse  mesh),  but  only  at  the  points  Vq-  (which  were  chosen  to  be  grid  points). 
For  first-order  difference  approximations  in  time,  this  is  equivalent  to  linear 
interpolation  in  time  between  n(>„-./„)  and  u(  r„*./„  -i-  A/).  Since  the  mesh 
ratio  AT,/(Ay, )‘  is  stable  (see  below),  we  expect  a  boundary  condition  obtained 
via  the  coarse  mesh  approximation  using  a  smaller  time  step  to  be  stable  as  well. 
Such  a  result  has  been  proved  in  other  contexts,  e.g.  [3]. 

As  has  already  been  indicated,  the  “rescaled  solution"  (/,  is  computed  using 
the  forward  Euler  scheme  in  the  variables  y,.  t,.  To  preserve  the  numerical 
scheme  and  to  maintain  accuracy,  it  is  important  to  use  the  same  discretization 
for  M,  as  for  m,  i.e.,  to  set  Ay,  =  A.v  and  At,  =  At.  This  requires  introducing 
new  grid  points:  those  u.sed  for  u.  .spaced  A.v  apart,  determine  m,(-.0)  only  on  a 
mesh  of  width  A  ‘  A.v.  To  achieve  Ay,  =  Ax.  A  '  -  1  new  points  must  be 
introduced  between  each  pair  of  existing  ones.  (If  for  example  A  =  then  3  new 
points  are  added  to  each  interval  to  refine  the  mesh  from  4  A.v  to  A.v.  The 
condition  that  A  ‘  be  a  smaller  integer  is  imposed  to  make  this  step  easy  to 
implement.)  Linear  interpolation  in  space  is  used  to  assign  a  value  to  m,(-.0)  at 
the  new  grid  points.  Notice  that  the  effect  of  this  procedure  in  the  original 
variable  x  is  to  refine  the  mesh  by  a  factor  of  A  '. 

The  structure  of  our  algorithm  should  now  be  clear:  after  rescaling,  the 
“original”  solution  u  and  the  “rescaled”  solution  h,  are  stepped  forward  inde¬ 
pendently,  each  on  its  own  grid.  A  single  time  step  of  u  corresponds  to  A  -  time 
steps  of  M,.  The  two  solutions  interact  primarily  in  that  u  is  u.sed  to  determine 
the  boundary  conditions  for  «,.  In  addition,  at  each  successive  time  step  for  ii 
the  coarse  grid  solution  is  modified  on  the  interval  that  was  re.scaled.  to  make  it 
agree  with  the  more  accurate  fine  grid  solution  When  ||«,(*,  ,VAt|)||^  first 
exceeds  M.  a  smaller  time  step  (equivalently,  linear  interpolation  in  time)  is  used 
to  find  a  value  t,*.  with  ( iV  -  1)  At,  ^  t,*  ^  (VAt,  such  that 

ll“i('- V)IL  =  A/. 
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On  the  interval  where  u,  >  aM  the  solution  is  rescaled  further,  yielding  m,.  and 
so  forth. 

The  (k  +  l)-st  rescaled  solution  is  introduced  when  the  “clock"  for 
reaches  a  value  such  that 

(2-13)  T;,*)l|^  =  A/. 

The  interval  in  y^-space  to  be  rescaled,  ( y*  ,  yy" ),  consists  precisely  of  the  grid 
points  where  aM: 

«A  (  >’a  -  V’a  .  V  )  <  «  ^  “a  {  -'  a  ■  T**  ) . 

(2.14) 

Ma  (  Vk  .  T**  )  ^  «A/  >  (  !•;  +  Ay^  .T*). 

The  next  rescaled  solution  ^ ,  is  related  to  by 

(2.15)  ^A  w( -'^A -  1' ’’’a  •  1 )  ~  ^  '**'a  (  ^-'a  •  1  • '’’a*  ^"■’"a  •  1 )  • 

its  “rescaled  space"  and  “rescaled  time"  variables  y^ .  ,  and  .  j  range  over 

(2.16)  X  '  'y^  ^  y* ,  i  ^  X  'y/  .  .  i  S  0. 

Its  initial  data  u^  ^  i(  vy,  i.O)  are  determined  by  rescaling  m^(  *.  ).  using  linear 

interpolation  to  define  it  on  a  refined  spatial  grid  of  mesh  width  A  v*  ,  j  =  A.v.  Its 
boundary  data  i(X  (vy* ,  /VAt^.,)  are  determined  from  using  the  coarse 
mesh  approximation,  and  it  is  stepped  forward  in  time  by  the  forward  Euler 
difference  scheme  with  At^  . ,  =  At.  Previously  rescaled  solutions  are  stepped 
forward  independently:  if  for  example  X  =  i,  then  is  stepped  forward  once 
every  16  time  steps  of  +  m*  i  once  every  256  time  steps  of  .  j,  etc. 
Whenever  a  fine  grid  solution  is  computed  at  a  point  in  space-time  where  a 
coarser  mesh  solution  is  also  defined,  the  value  of  the  coarse  mesh  solution  is 
updated  to  agree  with  the  fine  grid  calculation.  When  a  lime  step  is  reached  such 
that  A^At)||^  >  M.  then  it  is  time  for  another  rescaling.  A  smaller  time 

step  is  used  to  find  such  that  ||M^  i( ’•  ■’"a*  i )llx  X/.  and  the  entire 

procedure  is  repeated. 

This  algorithm  cannot  be  continued  indefinitely  without  losing  accuracy.  We 
use  the  symmetry  of  the  computed  solution  as  an  indicator  of  the  accuracy  of  the 
calculation.  With  symmetric  initial  data  and  perfect  arithmetic,  would  remain 
symmetric  for  every  k.  However,  the  roundoft  error  is  not  symmetric.  In  a  ty  pical 
calculation  using  X  =  T  the  amplitude  of  the  asymmetry  approximately  doubles 
from  one  rescaling  to  the  next.  On  a  Cray  XMP  using  double  precision  arith¬ 
metic,  machine  epsilon  is  about  10  and  by  the  87-th  iteration  only  two  or 
three  digits  of  symmetry  remain.  Generally  we  stop  the  calculation  after  80 
iterations;  with  p  =  5  and  X  =  ‘.the  corresponding  amplitude  of  u  is  on  the 
order  of  (X  =  lO’’.  The  u.se  of  multiple  grids,  rescaling  and  retming 
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only  where  the  solution  is  large,  is  crucial  to  the  success  of  this  calculation:  to 
achieve  similar  results  without  mesh  refinement  on  a  single,  uniform  grid  would 
require  10“^  mesh  points!  Similar  accuracy  could  be  achieved  using  a  single 
nonuniform  grid  (c.f.  [24],  [25])  but  our  method  has  the  advantage  of  choosing 
the  proper  distribution  of  grid  points  automatically. 

Certain  qualitative  features  of  the  .solution  are  strikingly  clear  from  even  a 
casual  examination  of  the  output  of  our  algorithm: 

(2.17a)  The  "rescaling  times”  t*  are  eventually  almost  independent  of  k 
(Figure  2). 

(2.17b)  If  i/*  ( • .  )  is  graphed  over  a  fixed  interval,  say  -  1  <  r  <  1 .  then  the 

graph  is  eventually  independent  of  k  (Figure.s  3.4). 

(2.17c)  The  width  of  the  interval  rescaled  at  the  A-th  iteration  beha\es  as  the 
.square  root  of  a  linear  function  (Figure  5). 

All  these  as.sertions  are  explained  in  Section  4.  It  turns  out  that  (2.17a)  is  a 
con.sequence  of  the  (proved)  result  (1.3).  while  (2.17b,c)  are  related  to  the 
conjectured  a.symptotics  (1.4). 

3.  Review  of  the  Theory 

Before  interpreting  the  algorithm,  we  review  some  of  the  results  and  conjec¬ 
tures  concerning  the  blow-up  of  u.  Our  discussion  is  restricted  for  simplicity  to 
the  case  at  hand:  positive,  symmetric,  radially  decreasing  solutions  of  (1.1 )  on  an 
interval  with  a  Dirichlet  boundary  condition.  It  should  be  noted,  however,  that 
parts  of  the  theory  have  been  carried  out  in  much  greater  generality,  including 
space  dimensions  n  >  1  and  non-radial  solutions. 

Let  T  be  the  blow-up  time  of  m,  in  other  words. 

||m(*. /)  !|.^  =  m(0, /)  -»  cc  as  /  T /'. 

It  is  convenient  to  introduce  “similarity  variables”,  a  change  of  both  dependent 
and  independent  variables  defined  by 


w(^.r)  =  (r-  z)'^"'’  ‘’m(.v./). 


(3.1) 


I  =  v/7^. 

V  =  -  log(  T  -  t). 
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One  computes  that  u  solves  (1.1)  if  and  only  if 

(3.2)  w,  -  Wjj  +  -  H'/’ =  0. 

Since  s  ->  oo  as  / 1  7',  the  change  of  variables  (3.2)  converts  any  question  about 
the  blow-up  of  u  into  one  about  the  large  time  a.symptotics  of  h  . 

Upper  and  lower  bounds  are  known  for  the  blow-up  rate  of  u  (see  [11 IV.  they 
imply  that 

(3.3)  0  <  C  ^  ■^)  IL  ^  ^ 

for  some  positive  C,  C\  independent  of  5.  Rewriting  the  equation  (3.2)  as 

w  —  —  ( p  vty  ) ,  -I-  - 7-  w'  —  H-  ^  =  0 

'  p  '  t '  £  p  -  I 

with  p(|)  =  exp{  -  \i^},  it  is  easy  to  see  that 

*  l\U  *  2(7^""  ' 

decreases  as  s  increases,  and  indeed  that 

^£[h1  s  -  jw}pdi. 

It  is  thus  natural  to  expect  that,  as  s  -►  00.  w(^.  s)  should  tend  to  a  stationary 
point  H-^(^)  of  the  functional  £.  It  turns  out  that  the  only  stationary  points  of  £ 
are  the*  trivial”  ones  =  0  and  =  ±(p  -  11  ".  Since  the  nonposi¬ 

tive  stationary  points  are  ruled  out  by  (3.3).  we  are  led  (heuristically)  to  the  con¬ 
clusion  that,  as  V  -►  to.  w($.  .r)  should  tend  to  the  constant  function  m^(0  = 
{p  —  1)" '/</’“ 0  These  ideas  can  be  made  rigorous,  and  they  lead  to 

Theorem  3.1  [14],  [16].  As  s  -*  00.  h’(|.  j)  -*  { p  -  ])  ''''£’  ".  uniformly 
on  the  set  |||  <  C  for  any  C  >  0. 

Transformed  into  the  original  variables  via  (3.1)  this  is  equivalent  to  (1.3). 


Now  consider  the  profile  of  w  as  a  function  of  |  for  fixed  .r  »  1.  Evidently  it 
is  nearly  flat,  w  =  ( p  -  1)'  **.  on  an  interval  about  the  origin  which  grows 

with  but  it  decays  to  0  at  the  endpoints  ^  =  ±e'^^.  corresponding  to  .v  =  ±  1. 
The  form  of  this  profile  can  be  guessed  by  supposing  that  m’(|.  s)  ~  f(i/g{s)) 
for  some  functions  f(r\)  and  gfs);  cf.  [14j.  [20j.  To  obtain  the  expected  qualita¬ 
tive  behavior  we  suppose  that 

(3.4a)  g(.5)T£«  as  r -*  00. 

and 


(3.4b)  /(O)  =  (p  -  1)  /(h) -^0  as 


850 


M.  BERGER  AND  R.  V.  KOHN 


moreover,  for  reasons  that  will  emerge  shortly,  we  also  want  to  assume  that 
(3.4c)  g'(s)/g{s)  -*0  as  s  -*  oc. 

i.e.,  that  g  has  subexponential  growth.  Substitution  of  the  ansatz  into  the 
equation  yields 

(3.5)  -  -  g  Hs)r{r,)  +  ~  0. 

As  i  -»  00,  the  first  two  terms  tend  to  0  by  virtue  of  (3.4a. c),  leading  to  this 
first-order  equation  for  /; 

(3.6)  - /'•{r,)  ==  0. 

The  general  solution  is 

(3.7)  /(’l)  =  ((iP  -  1)  +  cV) 

in  which  c  >  0  is  an  arbitrary  constant  of  integration.  Notice  that  /  satisfies 
(3.4b)  regardless  of  the  choice  of  c.  Absorbing  the  constant  of  integration  into  the 
(unknown)  g  leads  to 

Conjecture  3,2.  w'(|.  5) -((/?- 1)  +  for  some  g{s) 

such  that  g  — >  00  and  g'/g  0  as  s  oc. 


This  is  consistent  with  Theorem  3.1  since  $"/g(s)~  — *  0  as  .v  oc  with 
III  <  C.  As  will  be  explained  in  Section  4,  it  is  borne  out  by  our  calculations,  and 
indeed  is  responsible  for  the  observed  stability  of  the  profile  of  m^.  (2.17b).  In 
terms  of  the  original  variables,  this  conjecture  says  that 


(3.8)  u(x,t)  ~  (r- t)'’  ‘l(p-l)  + 


(r-t)gMliog(r-/)l) 


So  far  we  have  guessed  the  profile  /,  but  not  the  spreading  rate  g.  A  formal 
procedure  for  determining  g  is  worked  out  in  jl3],  (14),  following  a  method 
introduced  in  [20].  The  idea  is  to  look  for  a  formal  expansion  in  pow'ers  of  .? 
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and  to  suppose  that  the  first  two  terms  of  (3.9)  interact  by  setting 

(3.10)  g(T)=c„.v''^ 

(The  expansion  (3.9)  would  appear  to  be  consistent  with  other  choices,  for 
example  gis)  =  c„s,  but  for  such  g  it  becomes  impossible  to  satisfy  the  con¬ 
sistency  condition  (3.16)  below.)  The  value  of  the  constant  c,,  is  determined  by  a 
consistency  condition  for  the  existence  of  /,,  as  we  now  explain.  Substitution  of 
the  ansatz  (3.9)  into  the  equation  (3.2)  gives  a  .sequence  of  equations  of  orders 
■s",  .r  .s  etc.  The  first  .says  that  /„  sati.sfies  (3.6);  we  may  choose  the  constant 
of  integration  to  be  1  in  (3.7),  .since  c„  is  as  yet  undetermined  in  (3.10): 

(3.11)  Jaiv)  =  ((;>-  1)  +  r)  '  "• 

The  order  s  '  equation  says  that 

(3.12)  -  -  <■„  %”  +  -H  -  Pf','  Vi  =  (>• 

Using  (3.11),  this  takes  the  form 

(3.13)  /,'  +  «(»?)/,=  ^^(r)) 
with 

(3.14)  a(T})  =  -  — ^  ^  ,  .  /?(7))  +  3(',  V„". 

V  [P  1  (  /7  -  1  )  -L  TJ" 

We  can  write  (3.13)  as 

e-  ^(‘’V,)'= 

where  A  is  any  indefinite  integral  of  n.  One  computes  that 

=  c;/i/7, 

hence  the  general  solution  of  (3.12)  is 

(3.15)  /i(i))  , 

•'1 

with  Cj  e  R  a  new  constant  of  integration.  The  rule  proposed  in  [20]  for  choosing 
c,|  is  to  require  that 


(3.16) 


/i  should  be  analytic  at  tj  =  0. 
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In  other  words,  the  coefficient  of  r  in  the  Taylor  expansion 

(3.17)  +  «,/  +  aj'  +  ■■■ 

should  vanish,  so  that  the  corresponding  term  /f Tj’logf r; )  is  absent  fri'in  (3.1.‘>). 
The  logic  behind  this  requirement  is  that  such  a  logarithmic  term  wcuild  be 
difl'erentiated  in  the  process  of  finding  /,,  /,.•  -  •.  and  would  eventuallv  lead  to  a 
term  that  is  infinite  at  rj  =  0.  Calculation  gives  that 

(3.18)  u,  =  0  m  (3.17)  ^ — r, 

(P  1 )" 

The  first  term  in  the  expansion  (.^.9)  is  now  entirely  determined.  Since  our 
numerical  results  are  insufficient  to  resolve  the  next  term,  there  is  no  need  to 
evaluate  the  integral  (3.15)  for  /,.  However,  we  shall  make  use  of  its  value  at 
7]  =  0.  From  (3.15)  and  (3.17). 

/i(0)  =  -  !./;r(o)u„  -  - 1/?(0). 
where  /?  is  given  by  (3.14).  C  alculation  gives 

-y/3(o)  =  -c„'/,;'(0)  =  ^(/’  -  n  ' 

We  are  thus  led  to  this  refinement  of  Conjecture  3.2; 

CoN;t:(  rt’RE  3.3  (13).  (14).  Asympiotualh  as  s  -•  x. 

(3.19a)  w( .V )  =  (  /^  -  1 )  '  "  1  +  '  '  +  C(  V  ' ). 

Moreover,  at  i  =  0  the  order  s  '  correelion  to  (3.19a)  is  yiren  hr 

(3.19h)  vv(0.  v)  =  (p  -  \)  '  "  1  +  2^,.'  '  +  C(v  '). 

Rewritten  in  terms  of  the  original  variables  m(  .v. /).  the  first  assertion  is 
precisely  (1.4). 

The  preceding  calculation  is  purely  formal:  we  know  of  no  proof  that  the 
expansion  (3.9)  can  be  continued  to  all  orders,  or  that  it  correctly  represents  the 
behavior  of  w.  However,  there  is  .some  theoretical  support  for  (3.10):  it  is  known 
(.see  [10])  that  m(|,  .v)  -*  0  as  .s  ->  x  with  |^|/(.v‘  ’)  -►  x.  This  implies  that  if 
C'onjecture  3.2  is  valid  for  some  function  g(.v),  then  g  ^  c.v'  ^  (A  similar  result  is 
proved  for  a  slightly  difFerent  boundary  value  problem  in  [14].)  As  we  shall 
explain  presently,  the  spreading  rate  g(.v)  is  linked  to  the  width  of  the  interval 
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rescaled  at  each  stage  of  our  algorithm;  moreover,  the  linear  growth  of  {  )~ 

asserted  in  (2.17c)  reflects  the  conjectured  behavior  g  =  c„.v' 

4.  Numerical  Results  and  Interpretation 

In  this  section  we  use  our  numerical  results  to  test  the  conjectures  described 
above.  The  interpretation  of  the  results  is  complicated  by  the  presence  of  two 
sources  of  error;  the  discretization  error  in  using  the  forward  Euler  scheme  on  a 
finite  grid,  and  the  asymptotic  error,  which  arises  since  the  conjectures  refer  only 
to  the  behavior  of  vt(|,  s)  as  s  -»  oc. 

All  the  runs  reported  here  use  ^(.v)  =  1  +  cosItta  )  as  the  initial  data,  p  =  ^ 
for  the  nonlinearity,  and  X  =  ',  M  =  Itfl .  a  =  1. 8/2/2  for  the  parameters  of 
the  algorithm.  These  values  are  typical  but  arbitrary;  other  choices  of  initial  data 
and  algorithm  parameters  lead  to  identical  conclusions  about  the  asymptotic 
character  of  the  blow-up.  The  values  of  \x  and  Ar  are  always  cho.sen  so  that 
Ar/(A.v)'  =  .25.  and  each  run  is  continued  for  eighty  rescaling  iterations.  Since 
the  algorithm  is  an  unusual  one.  combining  rescaling  and  grid  refinement,  we  do 
a  convergence  study  for  Aa  =  .02,  .01,  and  .005  corresponding,  respectively,  to 
100.  200.  and  400  points  in  the  initial  grid  for  m(a.O). 

In  order  to  relate  our  calculations  to  the  conjectures  it  is  necessary  to  connect 
the  computed  “rescaled  solutions”  with  the  "real"  solution  u(\.t) 

and  the  "solution  in  similarity  variables”  m(^.  .v).  This  may  at  tirst  glance  appear 
impossible,  since  is  related  to  m*  .  ,  (and  hence  ultimately  to  u)  by  the  scaling 

(2.15) .  which  takes  place  at  the  implicitly  defined  lime  t/.  while  u  is  related  to  u 
by  the  change  of  variables  (3.1).  which  involves  the  unknown  blow-up  time  T.  In 
fact,  however,  it  is  possible:  the  missing  link  is  Theorem  3.1.  which  relates 
(7  —  / )  to  the  magnitude  of  u.  asymptotically  as  i  T. 

The  tirst  task  is  to  express  t^)  in  terms  i>f  u(.\.  i).  the  .solution  of  (1.1). 

If  is  the  “real”  time  at  which  the  rescaling  from  u^  to  ii^  . ,  takes  place,  then 

(2.1 5)  gives 

(4.1)  r*  =  T,* +  AV  +  ••• 

where  t*  is  the  .scaled  time  at  which  u,  is  rescaled  to  create  u^.  ,.  Iteratum  of 

(2.15)  also  gives  a  formula  for  the  computed  rescaled  solution  u*  just  befc^re  the 
next  re.scaling; 

(4.2)  «*(.,*.  V) 

where  }\  is  the  spatial  variable  of  «*.  In  particular,  at  lime  the  amplitude  of  ii 
has  increased  by  a  factor  of  \ 

(4.3)  =  <v(0./J  =  A  "A/, 

and  si)  the  blow-up  time  is  T  =  lim^ 
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Our  qualitative  observation  (2.17a)  is  concerned  with  the  number  of  time 
steps  taken  on  the  grid  for  before  the  creation  of  .  i-  Let  us  call  this  number 
A^^(IOO)  when  the  initial  grid  has  100  points,  and  so  forth.  Recall  that 
represents  the  “scaled  time"  for  (see  (2.15)),  and  that  is  created  at 

=  r*.  Since  we  use  the  same  time  step  =  A/  for  every  k.  would  be 
(Ar)  in  the  absence  of  discretization  error.  To  understand  the  behavior 
as  k  00,  we  write  (1.3)  in  the  form 

(4.4)  {T- +  o(]). 

where  o(l)  represents  a  term  that  tends  to  0  as  A  -►  oo.  When  combined  with 
(4.3)  this  gives  an  asymptotic  formula  for  7  -  in  terms  of  k: 

(4.5)  {T- +  o{]). 

The  behavior  of  t/  as  A  x  is  determined  by  (4.1)  and  (4.5): 

Ti*  =  x  -  >k-i) 

=  A-^*((r- /,.,)- (7- rj) 

=  1)  +  0(1). 

so  that 

(4.6)  limr**  =  M'  '-  !)• 

A;  -•  .V  P  ^ 

Thus  the  number  of  time  steps  should  become  asymptotically  independent 
of  A.  This  is  precisely  what  happens:  Figure  2  shows  A/  as  a  function  of  A.  for 
computations  using  100,  200  and  400  grid  points  initially,  and  Table  1  gives  the 
values  of  Nf.  at  selected  values  of  A. 

The  preceding  argument  was  based  on  the  rigorous  result  (1.3).  More  detailed 
asymptotics  for  can  be  obtained  from  the  conjectured  behavior  (3.19b),  which 
gives  a  l/s  correction  to  the  amplitude  of  m(0,  .v)  as  .v  -»  oc.  To  this  end,  let 
=  -log(7  -  /^),  and  note  from  (4.5)  that 

(4.7)  s^  =  2tlogA|A  +  [log(p  -  1)  +  {p  -  l)log  M]  +  o(l). 

From  (3.19b)  we  have 

w(0.  .?*)'’“'  =  +  0{s,  -). 


e 


?e 


25 


30 


35 


The  computed  time  ,Vj  At  until  the  amplitude  ofuj  reaches  threshold  M.  plotted  ag 

initial  grid  has  (a)  100  points,  (b)  200  points,  (c)  4(X)  points  Curse  (d)  is  the  pred 

At  from  (4.9) 


k 

^*(100) 

N*(200) 

\(400) 

20 

120.98 

478.82 

1910.21 

30 

120.20 

475.73 

1897.88 

40 

119.84 

474.28 

1892.09 

50 

119.63 

473.45 

1888.76 

60 

119.50 

472.91 

1886.60 

70 

119.40 

472.53 

1885.09 

80 

119.34 

472.25 

1883.98 
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Using  (4.7)  and  the  definition  of  w,  (3.1).  we  obtain 


(7-- 


/>-!  M 

4/)|log  X|  / 


+  0(k  -). 


strengthening  (4.4).  This  leads  to 


(7  -  / 


P  -  1 


in  place  of  (4.5),  .whence 


(4.8) 


=  * 


P  -  1 


(A-2-1)  1  + 


p-\  1 


4p|log  A| 


^  +  0(k  =). 


We  are  thus  led,  by  neglecting  the  error  term  in  (4.8).  to  an  asymptotic 
formula  for  N^\ 

(4.9)  Nr-’  -  (A,)-  •  -  l)(l  +  3^  •  I 


The  graph  of  •  At  is  the  lowest  line  in  Figure  2;  it  deviates  from  the 

values  computed  using  the  finest  mesh  by  only  about  .05%.  As  another  test  of 
(4.9),  we  note  from  Table  1  that  when  400  initial  grid  points  are  used,  the  number 
of  time  steps  per  iteration  changes  by  only  2.62  between  A  =  60  and  k  =  80.  For 
our  choices  of  A.  p,  and  M,  and  with  At  =  .25(Aa:)’  =  6.25  •  10  (4.9)  gives 

^^ymp  _  =  2.25, 


in  approximate  agreement  with  the  computation. 

The  observed  values  of  A*  can  also  be  used  to  test  the  convergence  of  the 
calculations  as  Ax  -*  0.  Focusing  on  k  =  80,  we  see  that  the  error  £  = 
1^80  -  .  A/  is 

£(400)  =  6.73  • 

£(200)  =  2.94  •  Ar^oo  =  11.76  •  At^y,. 

£(100)  =  2.01  •  A/,00  =  32.16  •  At^^. 

reflecting  the  first-order  convergence  rate  in  time  of  the  diflTerence  scheme. 

Our  second  qualitative  observation  is  concerned  with  the  computed  rescaled 
solution  just  prior  to  the  next  rescaling.  «*(>’*,  t**).  It  is  defined  for  -A  ‘v*  ,  < 
\\  <  A  ‘y**_ ,,  an  interval  that  grows  with  k.  To  compare  the  profiles  for  different 
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values  of  k.  it  is  therefore  convenient  to  rescale  each  )  so  as  to  make  it 

be  defined  on  a  fixed  interval  -1  <  r  <  1.  (Beware:  this  last  rescaling  afreets 
space  alone,  not  the  values  of  u^:  it  has  nothing  to  do  with  the  rescalings  that  are 
done  in  the  course  of  the  algorithm.)  Since  by  .symmetry  j  =  r**  i.  we  are  led 
to  consider  the  “  rescaled  profile” 

(4.10)  2  -  l.T**).  -1<2<1. 

Our  observation  (2.17b)  asserts  that  this  function  is  asymptotically  independent 
of  k. 

In  fact,  as  we  shall  now  explain,  the  form  of  this  rescaled  profile  can  be 
predicted  from  our  cruder  Conjecture  3.2  alone.  Indeed,  suppose  that 

(4.11)  h(|,s)  =/(4/g(5))  +  o(l) 

with  /  and  g  as  in  (3.4a-c).  To  relate  m*(  17,  )  and  w.  we  combine  (4.2)  and 

(3.1): 


(4.12)  M,(  y^.V)  ‘’(r- tj  ' ‘’M(XSy(r 


1 . ; 


.  v. 


with  5^  =  -log(r  -  1^).  as  above.  Since  w  is  a  uniformly  continuous  function  of 
its  arguments  (see  for  example  [16]).  (4.5)  and  (4.12)  yield 

(4.13)  u,(y*.  V)  =  (p  -  I)''*'’  ”  -y, .  v  J  +  D. 

Substitution  of  (4.11)  into  (4.13)  yields  the  corresponding  prediction  for  (4.10). 

=  (P-  1)‘/*'’  '>''2X  V  ,/8(^a))  +  <>(\). 

The  asymptotic  behavior  is  evidently  governed  by  that  of  the  ratio  \\  i/g(.?^). 
This,  too  can  be  evaluated  by  substituting  the  ansatz  (4.11)  into  the  known 
relation  between  and  w.  (4.13).  In  fact,  ignoring  errors  of  computation  for  the 
moment,  we  set  Ajc  =  0  in  the  definition  of  ly'  j.  (2.14).  to  get 

(4-14)  i(>'*-i- =  «*- i(.V*  i-T**,)=aM. 

Therefore,  from  (4.13)  and  (4.14),  we  have 

«A/=  (p  -  1)‘/''’-"a//(v/^^A/<'’  ,/g(.^  ,))  +  o(l). 

so  that  .V*  |/g(.s*  1)  tends  as  X  -»  00  to  the  (positive)  root  f  of 
(4.15)  a  =  (p-l)'^''’  "  -f). 
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The  observed  ratio  i/g(-f*)  has  the  same  asymptotic  value,  since  bv  Taylor's 
theorem 


gih) 


1  ^ 


max 

I  SoS  'i 


g'{o) 

,) 


i); 


using  this,  the  known  asymptotics  of  (4.7),  and  the  growth  hypothesis  on  g. 
(3.4c),  one  easily  shows  that 


lim 

k  OC 


gi^k  i) 


=  1. 


We  conclude  that  if  m>  has  the  proposed  form  (4.11).  then  the  rescaled  profile  of 
M*,  (4.10),  is  asymptotically 


(4.16)  r  ^  -  l)‘^''’  ‘’A//(v//>  -  1  M'/-  "  ’A 

Notice  that  the  predicted  profile  (4.16)  depends  only  on  /.  not  on  the  “spreading 
rate”  g(i):  the  effect  of  g  has  been  washed  out  by  the  final  re.scaling  of  w^. 

Conjecture  3.2  asserts  not  only  that  w  has  the  asymptotic  behavior  f(i/g(s)) 
but  also  the  form  of  /; 

/(i))  =  ((/^-i)  +  r)  '  ^ 


The  root  f  of  (4.15)  is  easily  computed  to  be 

(4.17)  f  -  1)‘  \ 

and  .substitution  into  (4.16)  yields 

(4.18)  m*(zA V)  ~  +  («'  '’-1)^  ■'■]  '  "■ 

The  right  side  of  (4.18)  is  virtually  indistinguishable  from  the  output  of  our 
algorithm  after  sufficiently  many  iterations,  see  Figures  3  and  4.  A  quantitative 
study  of  the  convergence  rate  is  difficult,  since  in  the  computation  the  set  where 
u  >  aM  is  enlarged  to  include  an  additional  grid  point  on  either  side  when 
rescaling  and  refinement  is  done.  However,  the  qualitative  behavior  can  be  seen 
by  evaluating  (4.18)  at  z  =  (.  and  comparing  it  with  the  rescaled  computed 
solution  at  z  =  i  after  k  =  gQ  iterations,  using  100.  200  and  400  points  in  the 
initial  spatial  grid: 

left  side  of  (4.18)  using  100  grid  points  =  1 .8081 . 
left  side  of  (4.18)  using  200  grid  points  =  1 .8064, 
left  side  of  (4.18)  using  400  grid  points  =  1 .8047. 


right  side  of  (4.18)  =  1 .8. 
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Figure  3.  The  compuled  profile  of  rescaled  as  in  (4  10)  a(  A  =  HO.  and  the  predicted  prolile 
{4.1S).  They  coincide  to  within  plotting  resolution. 


Recall  also  that  the  right-hand  side  of  (4.18)  does  not  include  a  ]/k  correction 
term,  but  is  the  asymptotic  limit  as  k  -*  yo. 

Our  third  qualitative  observation  is  that  (,v/)‘  grows  linearly  with  k.  This  is 
linked  to  the  “spreading  rate”  g  in  (4.11),  since  the  discussion  leading  to  (4.16) 
shows  that 

(4.19)  .v;  =g(5j(f +  o(l)). 

where  f  is  the  root  of  (4.15).  Conjecture  3.3  asserts  that 

gih)^  <0  =  7  ■ 

ip  - 


combined  with  the  asymptotics  of  5*.  (4.7),  and  the  value  of  (4.17).  this  yields 

(4.20)  =  ,;|logA|A/'  '’(a*  1)  +  r;(A). 

(f’  -  1) 
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'  0;/  i.  .  .  .  I  ,  .  , 

‘  .Z  ■  0 ,  0  Z  '  z 

f'igurc  4,  The  computed  profiles,  of  h*  for  .selected  values  of  A.  each  rescaled  as  in  (4  10)  As  A 
increases  they  converge  to  the  protile  shown  m  Figure  }■. 


As  in  the  last  paragraph,  we  expect  the  right  side  of  (4.20)  to  be  an  expansion  in 
powers  of  A.  In  particular,  if  (3.19a)  is  used  in  place  of  (4.1 1 )  in  the  derivation  of 
(4.19),  then  o(l)  becomes  0(1/A)  and  (4.20)  becomes 

(4.21)  (v;)'  =  A - ^L|logA|W*  '’(a'  '’-!)  +  0(1). 

(Z’  -  I) 

Thus,  Conjecture  3.3  predicts  not  only  that  (,v\*)’  is  asymptotically  linear,  hut 
also  the  value  of  the  slope. 

For  making  a  quantitative  comparison  between  our  numerical  results  and  the 
conjectured  behavior,  there  is  a  slight  advantage  to  replacing  i*'  in  the  above 
arguments  by  the  point  y*  defined  by 

(4.22)  u,{yt.T*)  =^0M. 


with  6  cho.sen  .so  that  a  <  <  1.  The  reason  is  that  the  numericallv  computed 

does  not  .satisfy  (4.14)  exactly,  since  it  is  required  to  be  a  grid  point,  see  (2.14). 
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Figure  5.  The  graph  of  (  ,v/)’  against  k.  for  various  values  of  9.  based  on  data  obtained  using  40<l 
points  in  the  initial  grid.  The  width  of  the  interval  rescaled  at  the  A  th  iteration  is  2  ri’.  with  o  as  in 
(2.2), 


However,  vf  need  not  be  a  grid  point,  so  it  can  be  chosen  to  solve  (4.22)  exactly 
(within  the  accuracy  of  the  calculation)  by  using  linear  interpolation  in  space. 
Replacing  a  by  0  in  the  arguments  that  led  to  (4.21)  we  see  that  (  tf)'  is 
expected  to  grow  as 


(4.23)  (yf)' =  y/c  +  0(1),  y  =  -,|logA|A/'  ''(6'  ^’-I). 

(p  -  1)' 


Thus  the  points  (/c,(  vf  )^),  with  0  held  fixed,  should  approach  a  line  of  slope 
y.  Our  values  for  p.  M,  and  X  give  y  =  .021{6  “  -  1).  Figure  5  shows  (  rf )'  as  a 
function  of  k  for  several  values  of  6.  using  400  points  in  the  initial  giid.  For  a 
quantitative  compari.son,  we  present  in  Table  2  the  slope  of  the  line  which 
best  fits  the  points  (X,(  vf)^)  in  the  least  squares  .sen.se  for  60  <  A  <  80.  and 
for  several  different  values  of  0.  The  results  agree  with  the  prediction  (4.23)  to 
within  3%. 
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Table  II.  Slope  of  line  through  the  points  (a,  (rf )’),  where  =  ff\f. 


OM 

computed 

slope 

predicted 

slope 

2.0 

.08327 

.08123 

2.25 

.04156 

.04054 

2.5 

.01772 

.01729 

2.75 

.00330 

.00322 

5.  Conclusions 

We  have  demonstrated  the  convergence  of  the  rescaling  algorithm,  and  used  it 
to  calculate  the  solution  of  u,  -  w,,  =  u-  until  its  magnitude  reaches  about  10''. 
The  computed  singularity  is  consistent  with  the  conjectured  form  (1.4).  derived 
by  means  of  a  formal  expansion. 

This  method  appears  to  be  .suitable  for  computing  singularities  that  arise  in 
the  solutions  of  other  equations  with  a  similar  scale  invariance.  Although  we 
obtain  satisfactory  results  using  the  forward  Euler  scheme,  which  is  just  hrst-order 
accurate  in  time,  it  may  be  better  for  some  applications  to  use  a  more  accurate 
discretization  of  the  partial  differential  equation.  A  natural  candidate  for  fur¬ 
ther  investigation  is  the  nonlinear  Schrodinger  equation  iu,  -  Aw  -  |m|'’  'i/  =  0 
on  a  ball  in  R"  with  radial  initial  data  and  a  Dirichlet  boundary  condition,  in  the 
critical  p  =  (n  ■¥  4)/«  or  supercritical  /?>(«  +  A)/n  cases.  Ihough  extensive 
calculations  have  already  been  done  (see  [24).  [25],  [27]).  we  think  that  our 
algorithm  may  achieve  a  greater  resolution  of  the  local  behavior  of  the  singularity 
than  those  done  to  date. 
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